/*
 * The MIT License
 *
 * Copyright 2013-2014 Florian Barras.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */
package jeo.math.linearalgebra;

import java.io.Serializable;

/**
 * Singular Value Decomposition.
 * <p>
 * For an m-by-n matrix A with m >= n, the singular value decomposition is an
 * m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and an n-by-n
 * orthogonal matrix V so that A = U*S*V'.
 * <p>
 * The singular values, sigma[k] = S[k][k], are ordered so that sigma[0] >=
 * sigma[1] >= ... >= sigma[n-1].
 * <p>
 * The singular value decompostion always exists, so the constructor will never
 * fail. The matrix condition number and the effective numerical rank can be
 * computed from this decomposition.
 * <p>
 * @author JAMA, http://math.nist.gov/javanumerics/jama/
 * @version 1.0.3
 */
public class SingularValueDecomposition
	implements Serializable
{
	////////////////////////////////////////////////////////////////////////////
	// ATTRIBUTE(S)
	////////////////////////////////////////////////////////////////////////////

	/**
	 * Generated serial version ID.
	 */
	private static final long serialVersionUID = 2256209583719000772L;
	/**
	 * Arrays for internal storage of {@code U} and {@code V}.
	 * <p>
	 * @serial internal storage of {@code U}
	 * @serial internal storage of {@code V}
	 */
	private final double[][] U, V;
	/**
	 * Array for internal storage of singular values.
	 * <p>
	 * @serial internal storage of singular values
	 */
	private final double[] s;
	/**
	 * Row and column dimensions.
	 * <p>
	 * @serial row dimension
	 * @serial column dimension
	 */
	private final int m, n;


	////////////////////////////////////////////////////////////////////////////
	// CONSTRUCTOR(S)
	////////////////////////////////////////////////////////////////////////////

	/**
	 * Constructs the singular value decomposition.
	 * <p>
	 * Note: structure to access {@code U}, {@code S} and {@code V}.
	 * <p>
	 * @param A a rectangular {@code Matrix}
	 */
	public SingularValueDecomposition(final Matrix A)
	{
		// Note: derived from LINPACK code
		// Initialize
		final double[][] components = A.getComponentsCopy();
		m = A.getRowDimension();
		n = A.getColumnDimension();
		/*
		 * Note: Apparently the failing cases are only a proper subset of (m<n),
		 * so let's not throw error. Correct fix to come later? if (m<n) { throw
		 * new IllegalArgumentException("Work only for m >= n"); }
		 */
		final int nu = Math.min(m, n);
		s = new double[Math.min(m + 1, n)];
		U = new double[m][nu];
		V = new double[n][n];
		final double[] e = new double[n];
		final double[] work = new double[m];
		final boolean wantu = true;
		final boolean wantv = true;
		// Reduce A to bidiagonal form, storing the diagonal elements
		// in s and the super-diagonal elements in e
		final int nct = Math.min(m - 1, n);
		final int nrt = Math.max(0, Math.min(n - 2, m));
		for (int k = 0; k < Math.max(nct, nrt); ++k)
		{
			if (k < nct)
			{
				// Compute the transformation for the k-th column and
				// place the k-th diagonal in s[k]
				// Compute 2-norm of k-th column without under/overflow
				s[k] = 0;
				for (int i = k; i < m; ++i)
				{
					s[k] = Norms.getEuclideanNorm(s[k], components[i][k]);
				}
				if (s[k] != 0.0)
				{
					if (components[k][k] < 0.0)
					{
						s[k] = -s[k];
					}
					for (int i = k; i < m; ++i)
					{
						components[i][k] /= s[k];
					}
					components[k][k] += 1.0;
				}
				s[k] = -s[k];
			}
			for (int j = k + 1; j < n; ++j)
			{
				if (k < nct & s[k] != 0.0)
				{
					// Apply the transformation
					double t = 0;
					for (int i = k; i < m; ++i)
					{
						t += components[i][k] * components[i][j];
					}
					t = -t / components[k][k];
					for (int i = k; i < m; ++i)
					{
						components[i][j] += t * components[i][k];
					}
				}
				// Place the k-th row of A into e for the
				// subsequent calculation of the row transformation
				e[j] = components[k][j];
			}
			if (wantu & k < nct)
			{
				// Place the transformation in U for subsequent back
				// multiplication
				for (int i = k; i < m; ++i)
				{
					U[i][k] = components[i][k];
				}
			}
			if (k < nrt)
			{
				// Compute the k-th row transformation and place the
				// k-th super-diagonal in e[k]
				// Compute 2-norm without under/overflow
				e[k] = 0;
				for (int i = k + 1; i < n; ++i)
				{
					e[k] = Norms.getEuclideanNorm(e[k], e[i]);
				}
				if (e[k] != 0.0)
				{
					if (e[k + 1] < 0.0)
					{
						e[k] = -e[k];
					}
					for (int i = k + 1; i < n; ++i)
					{
						e[i] /= e[k];
					}
					e[k + 1] += 1.0;
				}
				e[k] = -e[k];
				if (k + 1 < m & e[k] != 0.0)
				{
					// Apply the transformation
					for (int i = k + 1; i < m; ++i)
					{
						work[i] = 0.0;
					}
					for (int j = k + 1; j < n; ++j)
					{
						for (int i = k + 1; i < m; ++i)
						{
							work[i] += e[j] * components[i][j];
						}
					}
					for (int j = k + 1; j < n; ++j)
					{
						final double t = -e[j] / e[k + 1];
						for (int i = k + 1; i < m; ++i)
						{
							components[i][j] += t * work[i];
						}
					}
				}
				if (wantv)
				{
					// Place the transformation in V for subsequent
					// back multiplication
					for (int i = k + 1; i < n; ++i)
					{
						V[i][k] = e[i];
					}
				}
			}
		}
		// Set up the final bidiagonal matrix or order p
		int p = Math.min(n, m + 1);
		if (nct < n)
		{
			s[nct] = components[nct][nct];
		}
		if (m < p)
		{
			s[p - 1] = 0.0;
		}
		if (nrt + 1 < p)
		{
			e[nrt] = components[nrt][p - 1];
		}
		e[p - 1] = 0.0;
		// If required, generate U
		if (wantu)
		{
			for (int j = nct; j < nu; ++j)
			{
				for (int i = 0; i < m; ++i)
				{
					U[i][j] = 0.0;
				}
				U[j][j] = 1.0;
			}
			for (int k = nct - 1; k >= 0; --k)
			{
				if (s[k] != 0.0)
				{
					for (int j = k + 1; j < nu; ++j)
					{
						double t = 0;
						for (int i = k; i < m; ++i)
						{
							t += U[i][k] * U[i][j];
						}
						t = -t / U[k][k];
						for (int i = k; i < m; ++i)
						{
							U[i][j] += t * U[i][k];
						}
					}
					for (int i = k; i < m; ++i)
					{
						U[i][k] = -U[i][k];
					}
					U[k][k] = 1.0 + U[k][k];
					for (int i = 0; i < k - 1; ++i)
					{
						U[i][k] = 0.0;
					}
				}
				else
				{
					for (int i = 0; i < m; ++i)
					{
						U[i][k] = 0.0;
					}
					U[k][k] = 1.0;
				}
			}
		}
		// If required, generate V
		if (wantv)
		{
			for (int k = n - 1; k >= 0; --k)
			{
				if (k < nrt & e[k] != 0.0)
				{
					for (int j = k + 1; j < nu; ++j)
					{
						double t = 0;
						for (int i = k + 1; i < n; ++i)
						{
							t += V[i][k] * V[i][j];
						}
						t = -t / V[k + 1][k];
						for (int i = k + 1; i < n; ++i)
						{
							V[i][j] += t * V[i][k];
						}
					}
				}
				for (int i = 0; i < n; ++i)
				{
					V[i][k] = 0.0;
				}
				V[k][k] = 1.0;
			}
		}
		// Main iteration loop for the singular values
		final int pp = p - 1;
		// int iter = 0;
		final double eps = Math.pow(2.0, -52.0);
		final double tiny = Math.pow(2.0, -966.0);
		while (p > 0)
		{
			int k, kase;
			// Note: here is where a test for too many iterations would go
			/*
			 * This section of the program inspects for negligible elements in
			 * the s and e arrays. On completion the variables kase and k are
			 * set as follows:
			 *
			 * kase = 1 if s(p) and e[k-1] are negligible and k<p
			 *
			 * kase = 2 if s(k) is negligible and k<p
			 *
			 * kase = 3 if e[k-1] is negligible, k<p, and s(k), ..., s(p) are
			 * not negligible (qr step)
			 *
			 * kase = 4 if e(p-1) is negligible (convergence)
			 */
			for (k = p - 2; k >= -1; --k)
			{
				if (k == -1)
				{
					break;
				}
				if (Math.abs(e[k]) <= tiny + eps * (Math.abs(s[k]) + Math.abs(s[k + 1])))
				{
					e[k] = 0.0;
					break;
				}
			}
			if (k == p - 2)
			{
				kase = 4;
			}
			else
			{
				int ks;
				for (ks = p - 1; ks >= k; --ks)
				{
					if (ks == k)
					{
						break;
					}
					final double t = (ks != p ? Math.abs(e[ks]) : 0.) + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.);
					if (Math.abs(s[ks]) <= tiny + eps * t)
					{
						s[ks] = 0.0;
						break;
					}
				}
				if (ks == k)
				{
					kase = 3;
				}
				else if (ks == p - 1)
				{
					kase = 1;
				}
				else
				{
					kase = 2;
					k = ks;
				}
			}
			++k;
			// Perform the task indicated by kase
			switch (kase)
			{
				// Deflate negligible s(p)
				case 1:
				{
					double f = e[p - 2];
					e[p - 2] = 0.0;
					for (int j = p - 2; j >= k; --j)
					{
						double t = Norms.getEuclideanNorm(s[j], f);
						final double cs = s[j] / t;
						final double sn = f / t;
						s[j] = t;
						if (j != k)
						{
							f = -sn * e[j - 1];
							e[j - 1] = cs * e[j - 1];
						}
						if (wantv)
						{
							for (int i = 0; i < n; ++i)
							{
								t = cs * V[i][j] + sn * V[i][p - 1];
								V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - 1];
								V[i][j] = t;
							}
						}
					}
				}
				break;
				// Split at negligible s(k)
				case 2:
				{
					double f = e[k - 1];
					e[k - 1] = 0.0;
					for (int j = k; j < p; ++j)
					{
						double t = Norms.getEuclideanNorm(s[j], f);
						final double cs = s[j] / t;
						final double sn = f / t;
						s[j] = t;
						f = -sn * e[j];
						e[j] = cs * e[j];
						if (wantu)
						{
							for (int i = 0; i < m; ++i)
							{
								t = cs * U[i][j] + sn * U[i][k - 1];
								U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - 1];
								U[i][j] = t;
							}
						}
					}
				}
				break;
				// Perform one qr step
				case 3:
				{
					// Calculate the shift
					final double scale = Math.max(Math.max(Math.max(Math.max(Math.abs(s[p - 1]), Math.abs(s[p - 2])), Math.abs(e[p - 2])), Math.abs(s[k])), Math.abs(e[k]));
					final double sp = s[p - 1] / scale;
					final double spm1 = s[p - 2] / scale;
					final double epm1 = e[p - 2] / scale;
					final double sk = s[k] / scale;
					final double ek = e[k] / scale;
					final double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
					final double c = sp * epm1 * (sp * epm1);
					double shift = 0.0;
					if (b != 0.0 | c != 0.0)
					{
						shift = Math.sqrt(b * b + c);
						if (b < 0.0)
						{
							shift = -shift;
						}
						shift = c / (b + shift);
					}
					double f = (sk + sp) * (sk - sp) + shift;
					double g = sk * ek;
					// Chase zeros
					for (int j = k; j < p - 1; ++j)
					{
						double t = Norms.getEuclideanNorm(f, g);
						double cs = f / t;
						double sn = g / t;
						if (j != k)
						{
							e[j - 1] = t;
						}
						f = cs * s[j] + sn * e[j];
						e[j] = cs * e[j] - sn * s[j];
						g = sn * s[j + 1];
						s[j + 1] = cs * s[j + 1];
						if (wantv)
						{
							for (int i = 0; i < n; ++i)
							{
								t = cs * V[i][j] + sn * V[i][j + 1];
								V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + 1];
								V[i][j] = t;
							}
						}
						t = Norms.getEuclideanNorm(f, g);
						cs = f / t;
						sn = g / t;
						s[j] = t;
						f = cs * e[j] + sn * s[j + 1];
						s[j + 1] = -sn * e[j] + cs * s[j + 1];
						g = sn * e[j + 1];
						e[j + 1] = cs * e[j + 1];
						if (wantu && j < m - 1)
						{
							for (int i = 0; i < m; ++i)
							{
								t = cs * U[i][j] + sn * U[i][j + 1];
								U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + 1];
								U[i][j] = t;
							}
						}
					}
					e[p - 2] = f;
					// iter += 1;
				}
				break;
				// Convergence
				case 4:
				{
					// Make the singular values positive
					if (s[k] <= 0.0)
					{
						s[k] = s[k] < 0.0 ? -s[k] : 0.0;
						if (wantv)
						{
							for (int i = 0; i <= pp; ++i)
							{
								V[i][k] = -V[i][k];
							}
						}
					}
					// Order the singular values
					while (k < pp)
					{
						if (s[k] >= s[k + 1])
						{
							break;
						}
						double t = s[k];
						s[k] = s[k + 1];
						s[k + 1] = t;
						if (wantv && k < n - 1)
						{
							for (int i = 0; i < n; ++i)
							{
								t = V[i][k + 1];
								V[i][k + 1] = V[i][k];
								V[i][k] = t;
							}
						}
						if (wantu && k < m - 1)
						{
							for (int i = 0; i < m; ++i)
							{
								t = U[i][k + 1];
								U[i][k + 1] = U[i][k];
								U[i][k] = t;
							}
						}
						++k;
					}
					// iter = 0;
					--p;
				}
				break;
			}
		}
	}


	////////////////////////////////////////////////////////////////////////////
	// GETTER(S)
	////////////////////////////////////////////////////////////////////////////

	/**
	 * Returns the left singular vectors {@code U}.
	 * <p>
	 * @return the left singular vectors {@code U}
	 */
	public Matrix getU()
	{
		return new Matrix(m, Math.min(m + 1, n), U);
	}

	/**
	 * Returns the right singular vectors {@code V}.
	 * <p>
	 * @return the right singular vectors {@code V}
	 */
	public Matrix getV()
	{
		return new Matrix(n, n, V);
	}

	/**
	 * Returns the one-dimensional array of singular values.
	 * <p>
	 * @return the diagonal of {@code S}
	 */
	public double[] getSingularValues()
	{
		return s;
	}

	/**
	 * Returns the diagonal matrix of singular values {@code S}.
	 * <p>
	 * @return the diagonal matrix of singular values {@code S}
	 */
	public Matrix getS()
	{
		final Matrix X = new Matrix(n, n);
		final double[][] S = X.getComponents();
		for (int i = 0; i < n; ++i)
		{
			for (int j = 0; j < n; ++j)
			{
				S[i][j] = 0.0;
			}
			S[i][i] = s[i];
		}
		return X;
	}

	/**
	 * Returns the two norms.
	 * <p>
	 * @return max({@code S})
	 */
	public double norm2()
	{
		return s[0];
	}

	/**
	 * Returns the two norms condition number.
	 * <p>
	 * @return max({@code S})/min({@code S})
	 */
	public double cond()
	{
		return s[0] / s[Math.min(m, n) - 1];
	}

	/**
	 * Returns the effective numerical matrix rank.
	 * <p>
	 * @return the number of nonnegligible singular values
	 */
	public int rank()
	{
		final double eps = Math.pow(2.0, -52.0);
		final double tol = Math.max(m, n) * s[0] * eps;
		int r = 0;
		for (int i = 0; i < s.length; ++i)
		{
			if (s[i] > tol)
			{
				++r;
			}
		}
		return r;
	}
}
